The roles of stiffness and excluded volume in DNA denaturation 
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The nature and the universal properties of DNA thermal denaturation are investigated by Monte 
Carlo simulations. For suitable lattice models we determine the exponent c describing the decay of 
the probability distribution of denaturated loops of length I, P ~ l~ c . If excluded volume effects 
are fully taken into account, c = 2.10(4) is consistent with a hrst order transition. The stiffness of 
the double stranded chain has the effect of sharpening the transition, if it is continuous, but not of 
changing its order and the value of the exponent c, which is also robust with respect to inclusion of 
specific base-pair sequence heterogeneities. 
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The melting of a DNA molecule from a double stranded 
phase to a denaturated state where the two strands are 
unbound has been the subject of various studies in the 
past GrS|. Experiments done with UV absorption by 
diluted DNA solutions show that denaturation occurs 
through a series of jumps in the absorbance spectrum 
as a function of temperature H. The jumps are inter- 
preted as sudden denaturations of large double stranded 
portions of the inhomogeneous chain and suggest a first 
order character for the transition. The existing mod- 
els, even if believed to capture the relevant features of 
the system, introduce drastic simplifications of the com- 
plex structure of the DNA molecule in favor of analytical 
tractability. Within these models, two different mech- 
anisms responsible for a first order transition were re- 
cently proposed. For the Peyrard-Bishop (PB) model ||, 
it has been argued that the stronger stiffness of double 
stranded compared to single stranded DNA, may lead to 
a first order denaturation ||. It has also been proposed 
that this stiffness difference, in combination with base 
sequence heterogeneity, should be responsible for the ob- 
served jumps in absorption [ I| . Other studies using the 
Poland-Sheraga (PS) [l| model led to the claim that ex- 
cluded volume effects, even in the absence of stiffness, 
induce first order melting 0. 

In this Letter we investigate further these issues within 
models representing as realistically as possible the rele- 
vant properties of DNA. Within such models all mecha- 
nisms mentioned above can operate and thus be tested 
simultaneously without resorting to uncontrollable ap- 
proximations. Near the transition DNA can be regarded 
as an alternating sequence of double helix segments along 
which base pairs are bound, and of denaturated loops, 
where the two strands are detached. By a scaling analy- 
sis of the cumulative probability distribution function of 
denaturated loop lengths we show clearly that excluded 
volume effects drive the transition first order, in the limit 
of infinitely long chains. We also find that double helix 
stiffness alone does not modify such asymptotic result in 
a range of values chosen consistently with those expected 
for real DNA. However, in the presence of strong stiffness, 



and for a model with second order denaturation in the in- 
finite chain limit, the critical region becomes very narrow 
and, due to a slow crossover, thermodynamic quantities 
like the specific heat may behave consistently with first 
order even for long finite chains. 

We model the DNA strands by two self-avoiding walks 
(SAW) of length N on cubic lattice, identified by the 
vectors f\(i) and f^i) (0 < i < N), joined at a common 
origin (Fi(0) = ^(O)) and with free endpoints. A gain 
of energy e (= 1 here) is associated to a bond between 
the two strands, which occurs in the model when two 
monomers with the same i overlap (fi{i) — Tii%)). The 
binding energy is taken to be the same all along the chain, 
i.e., to start with, we neglect the heterogeneity of base 
pair interactions of specific sequences. At sufficiently low 
T the most probable configurations are fully bound, i.e. 
T\(i) = r*2(i) for all i. Upon increasing the temperature 
the two strands are expected to unbind at some T = T c . 
The transition is driven by the formation of denaturated 
loops, whose length can be measured by the number of 
unbound monomers, /, of the corresponding strand seg- 
ments. 




FIG. 1. Schematic representations of denaturated loops in 
model I and II. 

In order to investigate in detail excluded volume and 
stiffness effects we consider two different models, which 
we refer to as model I and II. While in both models 
bounded segments are self-avoiding and cannot overlap 
any other part of the chain, in I the excluded volume 
condition is partly relaxed (see Fig. Q): the two strands 
forming a denaturated loop are self-, but not mutually 
avoiding and overlap at non-complementary sites (i.e. 

= ^2 (fc) with j ^ k) are allowed. Such overlaps 
do not contribute to the energy. The stiffness of the dou- 
ble stranded chain is incorporated in a second parameter 
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Eb, which is the energy gain if two consecutive segments 
in the bound chain are aligned (i.e. when r\(i) = f^(i) for 
* = 3-h 3, i + 1 and n(j + 1) - n(j) = ri{i)-ri{j-l)). 
Of course, both models neglect the precise helix struc- 
ture. Model II, having excluded volume effects fully built 
in, is a more realistic representation of DNA. A recent 
study, in absence of stiffness || , claimed for it some evi- 
dence of first order transition. 

For models I and II we sampled by a multiple Markov 
chain Monte Carlo (MC) jlO method the loop length 
probability distribution function (pdf) at several differ- 
ent temperatures. In all cases the pdf has an exponential 
decay at low T (P(l, A) ~ exp(— I/Iq)) where bound seg- 
ments are instead broadly distributed in length. From 
a given T c upwards one observes instead a scaling form 
P(l,N) ~ l~ c f(l/N), where the exponent c seems to be 
approximately independent of T and / is a scaling func- 
tion. The length of bound double segments is narrow 
distributed in this case. The general picture emerging 
from our numerical results has an immediate thermody- 
namical interpretation. Quantities usually analyzed at 
polymer conformational transitions are the energy per 
monomer or the specific heat, which for a chain of length 
A and T ~ T r scales as HI: 



C(A, T) ~ N^^h [(T - T C )JV*] 



(1) 



where h is a scaling function and <j> is the crossover expo- 
nent. For large A one has C max (A) = maxTC(A, T) ~ 
A 2 *" 1 , from which cj> can be deduced. The density of 
binding contacts along the strands, proportional to the 
energy, should scale as A"^ 1 at T = T c . Clearly <f> < 1, 
and only <p = 1 implies a first order discontinuity of the 
density. Given the scaling form of P(l,N) and the fact 
that bound segment lengths are finite on average at T c , 
the same density should also scale as the reciprocal of 
the average loop length (I) = Y^i IP(1,N). Now, for 
1 < c < 2, (l) ~ A 2 ~ c , so that <j> = c - 1 < 1, and 
the transition is second order. If instead c > 2, (I) and 
the energy density remain finite at T c for A — > oo and 
the transition is first order (</> = 1). Analyzing P of- 
fers both fundamental and practical advantages over the 
use of Eq. (|l|). Indeed, due to finite size corrections to 
scaling, a reliable estimate of <fi requires good determina- 
tions of C around T c for sufficiently long chains |Q . On 
the contrary, the scaling behavior of P sets in already 
for relatively short chains (A rj 100), which allow us to 
estimate c reliably. The robustness of the estimate with 
respect to the variation of chain lengths assures that re- 
sults are little affected by finite size corrections. 

We consider first £& = 0. Figure |^ shows log-log plots 
of P as a function of I for A = 50, 100 and 180 for model 
I, at T = 0.85 w T c (we estimate T c = 0.85(2)). Af- 
ter an initial transient the data follow nicely a straight 
line in the plots, except when / Rj A, where of course 
P drops. It is interesting to note that the power-law 
regime sets in already for relatively short loops. A lin- 
ear fit of the data for A = 180 gives c = 1.73(4), indi- 
cating a continuous denaturation for model I. By fitting 



C max as C max (A) = AN 2 *" 1 (inset of Fig. |), we ob- 
tain <f> = 0.77(2). Thus again <j> < 1 indicates a contin- 
uous transition; moreover, the above proposed relation 
6 = c — 1 is well satisfied. 
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FIG. 2. Plot of P vs. I on a log-log scale for model I at 
T = T c and various chain lengths N. Inset: C max as a function 
of N and relative fit ~ jV^" 1 (see Eq. (@)). 

Within a PS model description each loop would be as- 
sumed to be totally uncorrelated with the rest of the 
chain, which at T > T c would just provide a grand- 
canonical critical reservoir of strand elements to form the 
loop. The value of c within this scheme can be calcu- 
lated easily. The critical grand-canonical pdf of a SAW 
of length I with fixed end-to-end distance R ei in three 
dimensions, is [01: 



p(l, R e 



F- 



-3u- 



(2) 



where g is a scaling function, 7 is the entropic, and v is 
the metric exponent. A loop in model I is made of two 
strands with common end-points and allowed to inter- 
sect each other. Thus the probability of a loop of length 
I (total perimeter 21) is: 



P(l) - / d 3 R eP 2 (l,R e )^l 2 ^ 



3u-2 



(3) 



from which follows c' IL ) = 2 + 3v — 27, where IL stands 
for isolated loop. Using the appropriate SAW exponents, 
i.e. 7 w 1.158, v Rj 0.588 Q, we find C ( IL ) Rj 1.448 p|. 
This result allows us to quantify the variation of the ex- 
ponent c due to excluded volume interactions between 
loops and segments which are taken into account in the 
MC simulations, namely Ac = c — c' IL ' Rj 0.3. 

We now turn to model II. Fig. || shows the loop proba- 
bility distribution at T Rj T c ( T c = 0.76(2)) for A = 50, 
100, 150 and 200. From a linear fit of the data we obtain 
c = 2.10(4) > 2. Excluded volume effects appear to be 
responsible for the discontinuous nature of denaturation, 
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in agreement with a recent extension of the PS model in 
which the interaction between loops and segments was 
included in an approximate way ffl. By rel yin g on field 
theoretical results for polymeric networks [fL4| , the au- 
thors of Ref. [Q estimated analytically the loop length 
pdf exponent to be c ~ 2.12, which is remarkably close 
to our numerical determination, suggesting that the sort 
of perturbative treatment of Ref. catches the essen- 
tial contribution to the correlations among different loops 
and segments. The inset of Fig. || shows a plot of P for 
T = 1.0, i.e. well above T c , and TV = 150; the data still 
show a power-law behavior with an exponent c = 2.09(5), 
i.e. consistent with the value found at T c . 
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FIG. 3. Log-log plots of P vs. I for model II at T = T c and 
for various TV. Inset: Same plots for T = 1.3T C . 



TABLE I. Summary of exponents found for the two models. 



Model 


C (IL) 


c(e b = 0) 


c(e b = 5) 


I 


1.448 


1.73(4) 


1.70(6) 


II 


1.762 


2.10(4) 


2.06(6) 



Table | summarizes the results obtained for the two 
models. In three dimensions the exponent of an isolated 
loop is c (IL) w 1.762 0. For both model I and II the ex- 
cluded volume interactions are responsible for a roughly 
identical increase of c with respect to its "bare" , isolated 
loop value c^ IL \ For model I this effect is not strong 
enough to cause a first order transition. 

Next we ask whether a sufficiently strong Eb > may 
induce first order denaturation. Recently such a pos- 
sibility has been discussed in the context of the PB 
model |,|. This model, like the PS one, takes into ac- 
count excluded volume effects very inadequately, since 
the two strands, while prevented from overlapping each 
other, are in fact not embedded in ordinary three dimen- 
sional space. Only a sort of longitudinal coordinate along 
the DNA backbone and the distance between strands en- 
ter the description. Within this context it was predicted 



that the stiffness of bound segments sharpens a continu- 
ous transition, making it look like first order for practical 
purposes Q] , or even drives it strictly first order || . 
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FIG. 4. Log-log plot of P vs. I for model I with Et = 5. 
Inset: scaling of C max as a function of TV". 

We focus first on model I, which has a continuous tran- 
sition for Eb — 0. Figure [| shows a plot of In P vs. In / 
for £b = 5. From a linear fit we obtain c = 1.70(6), in 
agreement with the value determined for Sb = 0. This 
indicates that a value of stiffness £& = 5 does not change 
the asymptotic character of the transition, which remains 
second order and in the same universality class as for 
Eb = 0. The inset of Fig. ^ shows a plot of C max as a func- 
tion of TV. This quantity scales linearly up to TV w 200, 
and deviates from linearity only for the longest lengths 
investigated. A linear scaling would imply (f> = 1, i.e. 
a first order transition. Apparently there is a very slow 
crossover in the specific heat, i.e. the deviation from 
Cmax (TV) ~ TV can be observed only for very long chains. 
Indeed, for e& = 10, which can also be a realistic value 
for DNA, the crossover must occur at TV values which 
are not accessible to our simulation. Notice that, on the 
contrary, already for rather short chains (TV < 100) the 
behavior of P indicates clearly that c ~ 1.7 < 2. So, 
while the investigation of C in the presence of realistic 
stiffness would reveal the true asymptotic nature of the 
transition only for extremely long chains, c seems to be 
very little affected by crossover. This discrepancy is due 
to the different way in which finite size effects and tem- 
perature uncertainties affect the scalings of P and C, as 
discussed above. 

For model II it is not possible to detect crossover be- 
havior in C max , since the transition is already first or- 
der ((f> — 1) at Eb — 0. For Eb = 5 and Eb = 10, 
which, as discussed below, are realistic choices, we es- 
timate c = 2.06(6) and c = 2.04(8), respectively. These 
values do not deviate, within error bars, from that ob- 
tained at Eb = 0. However, a very slight systematic de- 
crease of c with increasing Eb cannot be excluded. This 
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decrease could be determined by a mild crossover phe- 
nomenon induced by Eb > 0. Double stranded DNA has 
a persistence length which is about 20 — 50 times longer 
than that of single stranded chains |lq] . According to our 
calculations thepersistence lengths of bounded segments 
for T « 0.8T C f§ are about 10 and 30 times those of the 
unbound segments for £& = 5 and Eb = 10, respectively. 
Therefore the above parameter choices should be consid- 
ered as rather realistic for DNA. A single loop attached 
to very stiff long segments (a typical configuration for our 
chains at low temperatures) should experience very little 
excluded volume interaction with them. However, close 
to denaturation more and more loops start forming and 
the whole chain becomes rather flexible as the bubbles 
carry no stiffness. Typical estimates in our simulations 
yield a persistence length of just 2 — 3 lattice steps close 
to T c , explaining why the stiffness has little effect on the 
critical loop pdf. For a correct physical interpretation of 
the models studied here one should assign to each lattice 
monomer about 10 base pairs (bp), thus we predict that 
the persistence length of the double stranded chain close 
to denaturation would be of the order of 20 — 30 bp. 

We finally considered the effect of heterogeneity of the 
binding energies of base pairs along the chain. As an ex- 
ample we embodied in model II the information concern- 
ing a specific DNA sequence, resorting to a strictly micro- 
scopic interpretation of lattice monomers as single bases. 
We took Ecg/^at = 1.5 and Sb = 5, as expected for 
real DNA. For two different sequences of length N = 150 
fl7f we estimated c = 2.10(8), indicating that c robustly 
maintains the value determined for the homogeneous ver- 
sion of the model. 

Another important consequence of our results concerns 
codes |1^| used to simulate melting curves of real DNA 
sequences. These codes contain various parameters, such 
as the stacking energy of neighboring base pairs, and use 
c in the partition function for denaturated loops. The 
typical choice is c w 1.7 ]l8| , which is the appropriate 
value for a single isolated loop, while a more consistent 
choice would be c ~ 2.1. It would be interesting to in- 
vestigate the consequence of this different value of c on 
the calculated melting curves. 

Summarizing, the DNA denaturation transition can be 
characterized in terms of the algebraic decay of the cu- 
mulative loop length pdf. The exponent c signals clearly 
the asymptotic nature of the transition and is in fact 
well defined already for relatively short chains. In all our 
calculations we observe an onset of power-law behavior 
for short loops (/ > 5). This suggests that, if direct mea- 
surements could be realized, the sampling of P would not 
need to include very long denaturated loops and rather 
microscopic probes could be adequate. Promising in this 
respect seem to be techniques based on fluorescent DNA 
probes (see e.g. Ref. Q). 

Excluded volume effects alone appear to be responsi- 
ble for the discontinuous nature of denaturation in the 
infinite chain limit. Our calculations provide in fact a 
final verification of a conjecture advanced many decades 



ago by Fisher |2|] and recently corroborated and made 
more precise by Kafri et al. j?J . The stiffness of the dou- 
ble stranded DNA is not responsible for the first order 
character, contrary to claims based on results for the PB 
model H , but possibly only produces very mild crossover 
effects on c. In realistic conditions these effects should 
be barely detectable even on relatively short chains. 

Financial supports by MURST-COFIN 01, INFM- 
PAIS 01 and European Network ERBFMRXCT980183 
are gratefully acknowledged. 



[1] 
[3 



[4 
[5 

[6 

[7 

[8 

[9 

[10 

[11 

[12 

[13] 



[14] 
[15] 
[16] 



[17] 

[18] 
[19] 



D. Poland and H. A. Sheraga, J. Chem. Phys. 45, 1456 
(1966); 45, 1464 (1966). 

M. E. Fisher, J. Chem. Phys. 45, 1469 (1966). 
M. Peyrard and A. R. Bishop, Phys. Rev. Lett. 62, 2755 
(1989); T. Dauxois, M. Peyrard and A. R. Bishop, Phys. 
Rev. E 47, R44 (1993). 

D. Cule and T. Hwa, Phys. Rev. Lett. 79, 2375 (1997). 
N. Theodorakopoulos, T. Dauxois, and M. Peyrard, 
Phys. Rev. Lett. 85, 6 (2000). 

M. S. Causo, B. Coluzzi, and P. Grassberger, Phys. Rev. 
E 62, 3958 (2000). 

Y. Kafri, D. Mukamel, and L. Peliti, Phys. Rev. Lett. 85, 
4988 (2000). 

T. Garel, C. Monthus, and H. Orland, Europhys. Lett. 
55, 132 (2001). 

R. M. Wartell and A. S. Benight, Phys. Rep. 85, 67 
(1985). 

M. C. Tesi, E. J. Janse van Rensburg , E. Orlandini, and 
S. G. Whittington, J. Stat. Phys. 29, 2451 (1996). 
C. Vanderzande, Lattice Models of Polymers (Cambridge 
University Press, Cambridge 1998). 

The difficulty of a reliable determination of (f> from spe- 
cific heat data has been encountered in Ref. {J, where 
the estimate (f> ~ 0.9, for model II, is not fully consistent 
with the claimed first order nature of the transition. 
Notice that using the random walk exponents 7 = 1, 
v = 1/2 in the relation c = 2 + 3f — 27 one obtains 
c = 3/2, which is the correct value of c for a loop 
where no excluded volume effects are considered. It is 
also interesting to point out that for a loop in model I, 
c (il) _ 1 <- 3^2. This goes against the expectation 
that the incorporation of excluded volume effects always 
leads to an increase of the value of c. 
B. Duplantier, Phys. Rev. Lett. 57, 941 (1986). 
M. T. Record et al, Annu. Rev. Biochem. 50, 997 (1981). 
We assume a typical denaturation temperature of T c w 
80° C, therefore the room temperature, where persistence 
lengths are typically measured, corresponds to w 0.8T C . 
We used a DNA sequence from coliphage <I > X1 74 taken 



from http://www.ncbi.nlm.nih.gov/Taxonomy/ 



R. D. Blake et al, Bioinformatics 15, 370 (1999). 

G. Bonnet, O. Krichevsky, and A. Libchaber, Proc. Natl. 

Acad. Sci. 95, 8602 (1998). 



4 



